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Abstract. In numerically solving the Helmholtz equation inside a connected 
plane domain with Dirichlet boundary conditions (the problem of the quan- 
tum billiard) one surprisingly faces enormous difficulties if the domain has a 
problematic geometry such as various nonconvex shapes. We have tested sev- 
eral general numerical methods in solving the quantum billiards. Following 
our previous paper (Li and Robnik 1995) where we analyzed the Boundary 
Integral Method (BIM), in the present paper we investigate systematically 
the so-called Plane Wave Decomposition Method (PWDM) introduced and 
advocated by Heller (1984, 1991). In contradistinction to BIM we find that 
in PWDM the classical chaos is definitely relevant for the numerical accuracy 
at fixed density of discretization on the boundary b (b = number of numer- 
ical nodes on the boundary within one de Broglie wavelength). This can be 
understood qualitatively and is illustrated for three one-parameter families 
of billiards, namely Robnik billiard, Bunimovich stadium and Sinai billiard. 
We present evidence that it is not only the ergodicity which matters, but 
also the Lyapunov exponents and Kolmogorov entropy. Although we have 
no quantitative theory we believe that this phenomenon is one manifestation 
of quantum chaos. 
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1 Introduction 



It is quite embarrassing to realize that in an attempt to numerically solve 
the Helmholtz equation 

V>(r) + ArV(r) = 0, (1) 

satisfied by the scalar solution ^(r) with eigenenergy E = k 2 inside a con- 
nected plane domain B with the Dirichlet boundary condition ip(r) = on 
the boundary dB, one can face enormous difficulties in cases of "problem- 
atic" geometries such as e.g. various nonconvex shapes. This is precisely 
the problem of solving and describing the quantum billiard B as a Hamil- 
tonian dynamical system, which is thus just the 2-dim Schrodinger problem 
for a free point particle moving inside the enclosure dB, described by the 
wavefunction ip(r) with the eigenenergy E = k 2 . The corresponding classical 
problem is the classical dynamics of a freely moving point particle obeying 
the law of specular reflection upon hitting the boundary dB. Quantum bil- 
liards and their correspondence to their classical counterparts, especially in 
the semiclassical level, are important model systems in studies of quantum 
chaos (Gutzwiller 1990, Giannoni et al 1991). There are several general meth- 
ods for a numerical solution of equation ([]]) such as the Boundary Integral 
Method (BIM), recently carefully and extensively studied by Li and Robnik 
(1995a), following the important paper by Berry and Wilkinson (1984) (see 
also Boasman 1994), and the Plane Wave Decomposition Method (PWDM), 
introduced and advocated by Heller (1984,1991), whose analysis — especially 
in the light of the relevance of classical chaos — is the subject of our present 
paper. Another quite general method is the conformal mapping diagonaliza- 
tion technique introduced by Robnik (1984) and further developed by Berry 
and Robnik (1986) and recently by Prosen and Robnik (1993,1994), which 
in principle works for any shape whereas in practice it is used for shapes 
for which the conformal mapping onto the unit diskP] is sufficiently simple 
(possiblly also analytic). These methods can face quite similar problems in 
cases of almost intractable geometries, but they are to some extent comple- 
mentary. For example, the conformal mapping diagonalization technique can 
provide a complete set of all eigenenergies up to some maximal value beyond 
which the calculations cannot be performed due to the lack of computer stor- 
age (RAM), which means that we cannot reach very high- lying eigenstates. 

3 or some other integrable geometries admitting a simple basis for the representation 
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(Our present record (Prosen and Robnik 1994) is about 35,000 for the size of 
the banded matrix that we diagonalize in double precision, yielding at least 
12,000 good levels with accuracy of at least one percent of the mean level 
spacing.) 

However, using PWDM it is possible to go higher in energy by orders of 
magnitude but then only a few selected states can be calculated with many 
intermediate states in the spectral stretch missing. Therefore the geometry 
of some interesting and representative high-lying states can be analyzed, but 
the sample is typically not sufficiently complete (many missing states) to 
perform statistical analysis. See e.g. our recent papers on this topics (Li and 
Robnik 1994a, 1995b, 1995c, 1995d). The reasons for a failure of one of these 
methods can be quite different. For example, in BIM the main difficulty 
stems from the existence of "exterior chords" in nonconvex geometries in its 
standard formulation, but the trouble might be overcome by an appropriate 
reformulation of the method adapted to the correct semiclassical behaviour. 
This is discussed in (Li and Robnik 1995a) where we also show that classical 
chaos is completely irrelevant for BIM. On the contrary, in PWDM we find 
that the classical chaos is relevant for numerical accuracy especially in the 
semiclassical limit of sufficiently small effective Planck constant h e ff reached 
at sufficiently high eigenenergies. This demonstration and its qualitative 
explanation is the main subject of our present paper. To give a specific 
example we should mention isospectral billiards discovered and proved by 
Gordon et al (1992) which have been investigated experimentally (Sridhar 
and Kudrolli 1994) and it is also our experience (Li and Robnik 1994b) that 
BIM fails completely in this case (namely, due to strong nonconvexities) 
whereas PWDM at b = 12 yields the accuracy of eigenenergies of about 
within a few percent of the mean level spacing, except for some very special 
eigenmodes for which surprisingly we find agreement within double precision 
(16 digits) and which are characterized by the fact that these eigenvalues 
agree with the analytic solutions for the triangles within single precision (8 
digits). So the fact that in this and similar cases the experimental precision 
(for some levels) exceeds the best possible numerical precision even when 
using the best available methods is embarrassing for a theoretician but also 
motivation for further work. 
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2 The plane wave decomposition method of 
Heller 



In this section we present our general exposition of PWDM following Li 
and Robnik (1994a). To solve the Schrodinger equation ([!]) for ip(r) with 
Dirichlet boundary condition ip(r) = on dB we use the Ansatz of the 
following superposition of plane waves (originally due to Heller (1984)) 

N 

^( r ) = a j cos(k xj x + k yj y + fa), (2) 
i=i 

where k x j = kcosOj, k y j = k sin 9j, k 2 = E, and we use the notation 
r = (x,y). N is the number of plane waves and (j)j are random phases, drawn 
from the interval [0, 2ir), assuming uniform distribution, and 9j = 2j%/N de- 
termining the direction angles of the wavevectors chosen equidistantly. The 
Ansatz @) solves the Schrodinger equation (|l|) in the interior of the billiard 
region B, so that we have only to satisfy the Dirichlet boundary condition. 
Taking the random phases, as we discovered, is equivalent to spreading the 
origins of plane waves all over the billiard region, and at the same time this 
results in reducing the CPU-time by almost a factor of ten. For a given k 
we put the wavefunction equal to zero at a finite number M of boundary 
points (primary nodes) and equal to 1 at an arbitrarily chosen interior point. 
Of course, M > N. This gives an inhomogeneous set of equations which 
can be solved by matrix inversion. Usually the matrix is very singular, thus 
the Singular Value Decomposition (SVD) method has been invoked (Heller 
1984, Press et al 1986). After obtaining the coefficients aj we calculate the 
wavefunctions at other boundary points (secondary nodes). We always have 
3 secondary nodes between a pair of primary nodes. The experience shows 
that further increase of the number of secondary nodes does not enhance the 
accuracy. The sum of the squares of the wavefunction at all the secondary 
nodes (Heller called this sum "tension") would be ideally zero if k 2 is an 
eigenvalue and if (0) is the corresponding exact solution of (|I|). In practice 
it is a positive number. Therefore the eigenvalue problem now is to find the 
minimum of the "tension". In our numerical procedure we have looked for 
the zeros of the first derivative of the tension; namely the derivative is avail- 
able analytically/explicitly from (H) once the amplitudes aj have been found. 
In this paper we make the choice M = N, since it proved to be sufficient for 
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calculating the lowest 100 states whose accuracy we analyze. (For high-lying 
states studied in (Li and Robnik 1994a) we have used M = 5N/3.) It must 
be pointed out that the wavefunctions obtained in this way are not (yet) 
normalized, due to the arbitrary choice of the interior point where the value 
of the wavefunction has been arbitrarily set equal to unity. We therefore 
explicitly normalize these wave functions before embarking to the analysis of 
their properties. 

The accuracy of this method of course depends on the number of plane waves 
(N) and on the number of the primary nodes (M), and we have a considerable 
freedom in choosing iV and M > N. In order to reach a sufficient accuracy 
the experience shows that we should take at least N = 3£/ XdeBroglie, and 
M = N, where £ is the perimeter of the billiard and XdeBroglie is the de 
Broglie wavelength = 2n/k. With this choice in the present context and for 
the lowest 100 states we reach the double precision accuracy (sixteen dig- 
its) for all levels of integrable systems like rectangular billiard (where the 
eigenenergies can be given trivially analytically) and the circular billiard, 
but also for Robnik billiard B\ for small A < 0.1. Introducing the den- 
sity of discretization b defined as the number of numerical nodes per one de 
Broglie wavelength on the boundary we thus write the number of plane waves 
N = bC / XdeBroglie = b2n£/k. 

The main problem of investigation in this paper is to study the dependence 
of the systematical numerical error AE (i.e. the error due to the finite dis- 
cretization) on the density of discretization b, and the dependence of AE on 
the geometry (billiard shape parameter) at fixed b. In order to perform a 
systematic analysis the errors should be measured in some natural units and 
in our case this is of course just the mean level spacing which according to 
the leading term of Weyl formula is equal to An /A where A is the area of 
the billiard B. From now on we shall always assume that AE of a particular 
energy level is in fact measured in such natural units. Of course one imme- 
diately realizes that the error AE fluctuates wildly from state to state (see 
figure 5) so that generally nothing can be predicted about it individually. 
Therefore the approach must be a statistical one and so we typically take 
an average of the errors AE over a suitable ensemble of states. Specifically, 
in all cases of this paper we have taken the average of the absolute values 
of AE over the lowest 100 states (of a given symmetry class) and denote 
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it by < |Ai?| >. It is important and should be mentioned that we have 
also checked the stationarity of such average value over consecutive spectral 
stretches of 100 states each, so that our procedure does make sense. See a 
discussion at the end of section 3. 

It turns out that the accuracy of energy levels depends nontrivially on b, 
unlike in BIM (where we find always a power law (Li and Robnik 1995a)), 
namely it typically shows broken power law. By this we mean that 
< |Ai?| >obeys a power law < \AE\ >= Ab~ a with very large a for suffi- 
ciently small b, b < b c , whereas for larger b > b c it obeys a rather flat power 
law with very small positive a (close to zero). Therefore in contradistinction 
to BIM it is difficult to explore the general dependence of < \AE\ >on b, if 
there is any such universality at all. However, in order to investigate the de- 
pendence of the accuracy on geometry and the implied dynamical properties 
of billiards, we have decided to fix the value of b and have chosen b = 12, 
and then we look at the dependence of < \AE\ >on the shape parameters of 
three one-parameter billiards, namely Robnik billiard, Bunimovich stadium 
and Sinai billiard. 

3 Numerical results 

As is well known in classically integrable quantum Hamiltonian systems in 
the semiclassical limit (of sufficiently small h) the eigenfunction can locally be 
described by a finite superposition of plane waves with the same wavenum- 
ber, in case of plane billiards it is k — VE. If the quantum system has 
ergodic classical dynamics then in the semiclassical limit locally the wave- 
function can be represented as a superposition of infinitely many plane waves 
with the same k and with the wavevectors being isotropically distributed on 
the circle of radius k. Moreover, the ergodicity suggests to assume random 
phases for the ensemble of plane waves which implies that to the lowest ap- 
proximation the wavefunction is a Gaussian random function. While this is 
a good starting approximation, originally due to Berry (1977) and recently 
verified by Aurich and Steiner (1993) and also by Li and Robnik (1994a), 
the phases are actually not random but correlated in a subtle way dictated 
by the classical dynamics especially along the short and the least unstable 
periodic orbits, which is the origin of the scar phenomenon (Heller 1984, 
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1991, Bogomolny 1988, Berry 1989, Robnik 1989, Li and Robnik 1995d). 
(For a discussion see (Robnik 1988, 1995).) Thus we can qualitatively very 
well understand that PWDM should work well or even brilliantly in cases 
of classically integrable billiards whereas in the ergodic systems we expect 
a severe degradation of the accuracy (at fixed b) simply because the finite 
number of plane waves cannot capture the correct (infinite) superposition 
of plane waves everywhere in the interior of the billiard. If the system is a 
generic system of a mixed type with regular and irregular regions coexisting 
in the classical phase space, a scenario described by the KAM theory, then 
the degradation of accuracy (at fixed b) with increasing fractional measure of 
the chaotic component (denoted by p-i) is certainly expected. However, p2 is 
not the only parameter which controls the accuracy (at fixed b), since, as we 
shall see, the dynamical properties like diffusion time, Lyapunov exponent 
and the Kolmogorov entropy play a role. It is the aim of the present paper to 
numerically explore this type of behaviour in three different billiard systems. 

The first billiard system is defined as the quadratic (complex) conformal map 
w = z + \z 2 from the unit disk \z\ < 1 from the z plane onto the w = (x, y) 
complex plane. The system has been introduced by Robnik (1983) and fur- 
ther studied by Hayli et al (1987), Frisk (1990), Bruus and Stone (1994) and 
Stone and Bruus (1993a,b) for various parameter values A. Since the billiard 
(usually called Robnik billiard) has analytic boundary it goes continuously 
from integrable case (circle, A = 0) through a KAM-like regime of small 
A < 1/4 with mixed classical dynamics, becomes nonconvex at A = 1/4 (the 
bounce map becomes discontinuous), where the Lazutkin caustics (invariant 
tori) are destroyed giving way to ergodicity. As shown by Robnik (1983) the 
classical dynamics at these values of A is predominantly chaotic (almost er- 
godic), although Hayli et al (1987) have shown that there are still some stable 
periodic orbits surrounded by very tiny stability islands up to A = 0.2791. 
At larger A we have reason and numerical evidence (Li and Robnik 1994c) 
to expect that the dynamics can be ergodic. It has been recently rigorously 
proven by Markarian (1993) that for A = 1/2 (cardioid billiard) the system 
is indeed ergodic, mixing and K. This was a further motivation to study the 
cardioid billiard classically, semiclassically and quantanlly by several groups 
e.g. by Backer et al (1994) and Backer (1995), Bruus and Whelan (1995). 
The billiard shape for A = 0.4 is shown (the upper half) in figure la. Since 
all states are either even or odd we can take into account these symmetry 
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properties explicitly. In fact we want to specialize to the odd eigenstates only. 
Therefore in order to a priori satisfy the Dirichlet boundary condition on the 
abscissa of figure la we specialize the general Ansatz (g) to the following form 

N 

i/)(t) = a j cos(k xj x + (pj) sm(k yj y), (3) 
i=i 

where all the quantities are precisely as in (Q) except that the N discretiza- 
tion (primary) nodes are equidistantly located only along the half of the full 
billiard boundary, so that b is exactly the same as in using the Ansatz fl2|) for 
the full billiard. 

In figure 2 we show the results for this billiard, namely we plot < \AE\ > versus 
A at fixed b = 12. Close to integrability (A < 0.1) we reach the accuracy 
within 14 to 15 digits which is almost the double precision on our machine 
(16 digits), in which all our calculations have been performed. As the value 
of A increases we observe a dramatic deterioration of the accuracy where 

< \AE\ > increases by many orders of magnitude, namely almost by 13 
decades, leveling off at < \AE\ > approximately equal to 10 -2 , which means 
that we have now the accuracy of only a few percent of the mean level spac- 
ing. This dramatic but quite smooth increase of < \AE\ >is certainly related 
to the emergence of classical chaos with increasing A, but definitely is not 
controlled merely by p%, because p2 reaches the value of 1 (almost ergodic- 
ity) already at A = 1/4 (Robnik 1983, Prosen and Robnik 1993), whereas 

< \AE\ >still varies considerably in the region A > 1/4. Thus it is obvi- 
ous that in the semiclassical picture also other classical dynamical properties 
(measures of the "hardness" of chaos) play an important role. Although we 
do not have a quantitative theory yet one should observe that according to 
(Robnik 1983) the Lyapunov exponent and Kolmogorov entropy (h) vary 
also quite smoothly with A, suggesting a speculation that there might be a 
relation between < \AE\ >and h. 

Another demonstration of the effectivity of PWDM and its accuracy is dis- 
played in table 1 where we show the numerical value of the scalar product of 
two consecutive normalized eigenstates, namely the ground state and the first 
excited state, denoted by Oi 2 , which ideally should be zero. As we see here 
too the accuracy decreases (by orders of magnitude) sharply but smoothly 
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with increasing shape parameter A. 



It is then interesting to similarly analyze an ergodic system such as sta- 
dium of Bunimovich shown in figure lb, where the shape parameter is a/R 
and we have looked at the results for < a/R < 10. In fact for our purposes 
we have chosen and fixed R — 1 in all cases. Since this billiard is known to 
be rigorously ergodic (and mixing and K) for any a > in this case P2 is 
exactly 1 and constant. We have calculated the lowest 100 energy levels of 
the odd-odd symmetry class. Therefore in this case the general Ansatz @) 
can be specialized as follows 

N 

^( r ) = a i sin-(k xj x) sm(k yj y). (4) 

Here again the discretization (primary) nodes are only on the outer boundary 
of the stadium with discretization density b = 12. From our plot in figure 
2 we see that in the integrable case of the circle (a = 0) we again reach the 
accuracy within at least 14 digits, but this brilliant accuracy at fixed b = 12 
deteriorates almost discontinuously upon increasing a and then < | A.E| >still 
increases by about two orders of magnitude when a goes from 0.1 to 10. It 
appears to us that classical chaos is definitely relevant for the accuracy of 
the method which might and should be explained by an appropriate theory 
in the semiclassical level. As an observation we should mention that the Kol- 
mogorov entropy increases sharply with a/R when a/R goes from to about 
1, where it reaches the maximum, and then decreases slowly (Benettin and 
Strelcyn 1978), whereas our < \AE\ >increases monotonically. Thus if there 
is a relationship between < |A.E| >and Kolmogorov entropy it certainly is 
not a simple one. 

We have tested also another system with hard chaos, namely the Sinai bil- 
liard sketched in figure lc (desymmetrized). The system is known to be 
ergodic, mixing and K. In calculating the 100 lowest energy levels of the 
desymmetrized Sinai billiard we used the same specialized Ansatz as in eq. 
@, thereby taking into account explicitly the Dirichlet boundary condition 
on the abscissa y = 0. In this case b is the density of discretization of the 
equidistant nodes along the rest of the perimeter. Similarly as in case of 
stadium we easily reach the double precision of 16 digits in the limiting in- 
tegrable case of zero radius R = 0, but this accuracy is almost instantly lost 



9 



by increasing R as seen in figure 3. The value of < \AE\ >levels off at about 
lCT 3 to lCT 2 for all R between 0.025 and 0.45. 

As a final technical point we comment on the stationarity of < \AE\ >as 
a function of energy, which has been confirmed for the Robnik billiard at 
A = 0.27 where the average value over consecutive spectral stretches over 
100 states has been found to be quite stable. Specifically, to illustrate this 
finding we plot in figure 5 the absolute values of the lowest 400 consecutive 
eigenstates where one can see that the average value over 100 consecutive 
states is quite stable indeed. This is shown in table 2 for four intevals of 100 
states each. 

4 Discussion and conclusions 

We believe that our present paper presents quite firm numerical (phenomeno- 
logical) evidence for the relevance of classical chaos for the effectiveness of 
PWDM as a quantal numerical method to solve a quantum billiard, which 
is manifested especially in the semiclassical limit and might and should be 
explained in terms of an appropriate semiclassical theory. Qualitatively the 
reasons for this phenomenon are explained in the first paragraph of section 3. 
The parameter p 2 (the fractional volume of the chaotic component (s)) def- 
initely plays an important role but is not the only aspect of classical chaos 
controlling the behavior of the error < \AE\ >at fixed discretization den- 
sity b. Namely even in rigorously ergodic systems where p 2 is 1 the error 

< \AE\ >might be controlled by the slow diffusion in the classical phase 
space (diffusive ergodic regime, soft chaos). If the classical diffusion time 
is much longer than the break time threat {tbreak = h/D, D = mean en- 
ergy level spacing) then the quantal states will be strongly localized in spite 
of the formal ergodicity (for a demonstration see (Li and Robnik 1995b)), 
and therefore they mimic certain amount of regularity enabling a better ac- 
curacy of PWDM, i.e. < \AE\ >is smaller than for completely extended 
chaotic high-lying eigenstates where according to our experience somehow 

< \AE\ >typically saturates at about a few percent of the mean level spac- 
ing, even if we brutally increase b beyond any reasonable limits. Indeed, as 
can be seen by comparison of figures 2, 3 and 4, in case of the stadium this 
saturation value of < \AE\ >is about 10" 4 , which is almost two orders of 
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magnitude smaller than in the Sinai billiard (figure 4) and cardioid billiard 
(figure 2). We think that this is due to the strong localization of eigenstates 
in the stadium, which is very well known to display an unusual abundance 
of scars (Heller 1984, Li and Robnik 1995d). 

Thus our present work is a motivation for a semiclassical theory to explain 
this aspect of quantum chaos which exhibits some algorithmical properties of 
PWDM in applying it to quantum billiards with variety of classical dynamics. 
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Tables 



Table 1. The test of the orthogonality of the eigenstates. The scalar product 
of two consecutive normalized wavefunctions Ou, namely the ground state 
and the first excited state, for Robnik billiard at different shape parameters. 



A 


Ol2 





2.0E-16 


0.1 


-5.0E-15 


0.2 


7.8E-10 


0.3 


4.8E-6 


0.4 


5.5E-4 


0.5 


1.5E-3 
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Table 2. The stationarity test of < \AE\ >for Robnik billiard at A = 0.27 
for the lowest 400 odd eigenstates. 



average stretch 


< \AE\ > 


1-100 


1.54E-7 


101-200 


2.21E-7 


201-300 


2.77E-7 


301-400 


2.03E-7 


1-400 


2.13E-7 
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Figure captions 



Figure 1 (a-c): The geometry of the boundary of the three desymmetrized 
billiards, the Robnik billiard (a), the Bunimovich stadium (b) and the Sinai 
billiard (c). 

Figure 2: The ensemble averaged (over 100 lowest odd eigenstates) absolute 
error (measured in units of mean level spacing) versus the billiard shape 
parameter A for Robnik billiard with fixed density of boundary discretization 
b = 12. The numerical points are denoted by the diamonds, which are joined 
by straight lines just to guide the eye. 

Figure 3: The ensemble averaged (over 100 lowest odd-odd eigenstates) 
absolute error (measured in units of mean level spacing) versus the billiard 
shape parameter a/R for Bunimovich stadium with fixed density of boundary 
discretization b = 12. The numerical points are denoted by the diamonds, 
which are joined by straight lines just to guide the eye. 

Figure 4: The ensemble averaged (over 100 lowest eigenstates) absolute 
error (measured in units of mean level spacing) versus the billiard parameter 
R (the radius of inner circle) for the desymmetrized Sinai billiard with fixed 
density of boundary discretization b = 12. The numerical points are denoted 
by the diamonds, which are joined by straight lines just to guide the eye. 

Figure 5: The absolute error of eigenstates versus eigenenergy for the low- 
est 400 odd eigenstates of Robnik billiard at A = 0.27. The averages over 
consecutive stretches (of 100 states each) are given in table 2, demonstrating 
that < \AE\ >is quite stationary. 
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Figure 4 
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Figure 5 
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